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We study the mechanical stiffening behavior in two-dimensional (2D) cross-linked networks of semiflexible 
biopolymer filaments under simple shear. Filamental constituents immersed in a fluid undergo thermally excited 
bending motions. Pulling out these undulations results in an increase in the axial stiffness. We analyze this 
stiffening behavior of 2D semiflexible filaments in detail: we first investigate the average, static force-extension 
relation by considering the initially present undulated configuration that is pulled straight under a tensile force, 
and compare this result with the average response in which undulation dynamics is allowed during pulling, as 
derived earlier by MacKintosh and coworkers. We will show that the resulting mechanical behavior is rather 
similar, but with the axial stiffness being a factor 2 to 4 larger in the dynamic model. Furthermore, we study the 
stretching contribution in case of extensible filaments and show that, for 2D filaments, the mechanical response 
is dominated by enthalpic stretching. Based on the single-filament mechanics, we develop a 2D analytical model 
describing the mechanical behavior of biopolymer networks under simple shear, adopting the affine deformation 
assumption. These results are compared with discrete, finite-element (FE) calculations of a network consisting 
of semiflexible filaments. The FE calculations show that local, nonaffine filament reorientations occur that 
induce a transition from a bending-dominated response at small strains to a stretching-dominated response at 
larger strains. Stiffening in biopolymer networks thus results from a combination of stiffening in individual 
filaments and changes in the network topography. 

PACS numbers: 87.16.Ka, 87.15.La, 82.35.Lr, 82.35.Pq 



I. INTRODUCTION 

The mechanics of eukaryotic cells is largely governed by the 
cytoskeleton, an interpenetrated network of biopolymer pro- 
tein filaments, spanning the region between the cell's nucleus 
and membrane. Its filamentous constituents are actin fila- 
ments, intermediate filaments and microtubules. It is impor- 
tant to gain more insight in the physical origin of the me- 
chanical behavior of the cytoskeleton in view of its role in 
fundamental biological processes as cell division, cell motil- 
ity and mechanotransduction. It is well known that network- 
like biological tissues are compliant and respond to defor- 
mation by exhibiting an increasing stiffness, i.e., the ratio 
of change of stress and change of strain. Experimental re- 
search on the viscoelastic strain stiffening behavior in such 
biopolymer networks started in the 1980's, e.g. by means of 
micropipette and microtwisting experiments^ and rheological 
experiments on in-vitro gels of cytoskeletal filaments (actin, 
vimentin, keratin) , ^i^'^i^'^ neuronal intermediate filaments^ 
and fibrini^ Through these and numerous other studies, it 
is well-estabUshed by now that not only the filamental con- 
stituents but also the type and dynamics of cross-linking pro- 
teins govern the mechanical response of such semiflexible 
biopolymer networks i^^i^'i'^i^^''"*''^ 

The mechanical behavior of semiflexible biopolymer net- 
works has been subject to theoretical investigation since the 
mid nineties , i i"^i i ^- 1 9,20^2 1 ^22,23 ]y[ost studies focussed on the 
small-strain deformations enabling analytical treatment. The 
mechanics of these networks under strain is determined by the 
mechanical properties of individual filamental constituents, 
changes in network topography under deformation, and fi- 
nally, the molecular properties and dynamics of cross-binding 
proteins. MacKintosh and co-workers developed a model 
to describe the elasticity of biopolymer networks, in which 



strain stiffening originates from longitudinal stiffening in in- 
dividual semiflexible filaments.!^ Filamental stiffening is at- 
tributed to entropic effects, originating from the dynamic in- 
teractions with the surrounding (cytoplasmic) fluidj '^'^^'^^i^^ 
Next, Storm et al. developed a continuum model by insert- 
ing this single-filament force-extension relation in a network 
of infinitely many, randomly oriented filamentsi^ Under an 
applied shear, the network is then assumed to deform in an 
affine manner, allowing for an analytical description of the 
overall nonlinear network response. 

Under the affine-deformation assumption, used in network 
models describing rubber elasticity^^ strain stiffening of the 
network directly results from stiffening in individual fila- 
ments. The validity of the affinity assumption adopted by 
Storm et al.^ was studied in detail by Head et alJ^^ and by 
Wilhelm and Freyi^ for straight filaments. Numerical calcula- 
tions in the small-strain limit show that network distortions are 
indeed affine for high density, highly cross-linked networks, 
but are nonaffine for low- and intermediate-density networks. 
In a recent publication we have reported on strain stiffening 
in two-dimensional (2D) cross-linked biopolymer networks 
comprising discrete filaments, analyzed by the finite-element 
method)^ With the aid of this novel approach, the network 
response was calculated up to large strains for various fila- 
ment densities, and we showed that the origin of stiffening 
lies in the network rather than in its constituents. In this dis- 
crete filament model strain stiffening results from nonaffine 
network reorientations that induce a transition from a (rela- 
tively soft) bending-dominated response at small strains to a 
(stiff) stretching dominated response at large strains. Next to 
initially straight filaments, we studied networks of initially un- 
dulated filaments as if instantaneous thermal undulations are 
frozen in, and found that these merely postpone the onset of 
the stiffening. 
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Both approaches have advantages and disadvantages. The 
affine-network (AN) model accounts for the dynamic inter- 
action with the surrounding fluid during pufling, but neglects 
any network rearrangements, while the discrete-network (DN) 
model accounts for the latter, but ignores the former. The aim 
of this article is to rigorously investigate the relative contri- 
bution of both effects on the overall stiffening by a detailed 
comparison of both approaches. 

This article is organized as follows. In Sect. II we will 
determine the undulated shape of 2D semiflexible filaments 
resulting from the interaction with the surrounding fluid and 
will investigate the distribution of slack (the relative displace- 
ment of the ends to reach a straight configuration). Subse- 
quently, in Sect. Ill we will study the mechanical behavior of 
extensible and inextensible semiflexible filaments under ten- 
sion. We will in detail compare the ensemble-averaged force- 
extension relation of individual, inextensible filaments with 
and without undulation dynamics taken into account (to be 
referred to as the dynamic and static response, respectively) 
and show that the average mechanical response is rather simi- 
lar. In addition, we will investigate the average static and dy- 
namic force-extension relations in case of extensible filaments 
and show that the mechanical response is largely determined 
by the enthalpic stretching contribution, as confirmed by dis- 
crete finite-element calculations. Finally, in Sect. IV we will 
develop an analytical model describing the affine, mechanical 
behavior of 2D biopolymer networks under simple shear. As 
input, the model uses the average force-extension relation of 
individual filaments determined in Sect. III. Results will be 
compared with the discrete network model and we will show 
that nonaffine local network reorientations, that are not taken 
into account in the analytical model, play a key role in the 
onset of stiffening, especially at lower densities. 



II. SEMIFLEXIBLE FILAMENTS: SHAPE AND 
THERMAL UNDULATIONS 



Mathematically, the persistence length follows from the av- 
erage cosine of 9{s) — 0(0) in a filament over time or, alter- 
natively, over all filaments in the ensemble at a specific time. 
For 2D filaments this average decays exponentially with arc 
length s as^ 



(cos[0(s)-0(O)]) =exp 



with the persistence length defined by 



s 

2L^ 



Lp = 



(1) 



(2) 



In Eq. (|2]), is Boltzmann's constant, T is the temperature, 
and Kf is the flexural rigidity. While biopolymers such as DNA 
are flexible polymers with Lp ^ Lc, we here focus on semi- 
flexible filaments (e.g. actin, vimentin, keratin and fibrin) that 
have a persistence length similar to their contour length. De- 
tails on the values of Lp and Kf used for the calculations in this 
article can be found in Appendix A. 

In general, the undulated shape of filaments can be ex- 
pressed by a superposition of Fourier modes for the tangent 
angle OisyM^ 



= W — ^fl"cos [q„s] 



(3) 



with q„ = nn/Lc- The amplitudes aj] can be calculated for a 
given shape 9 (s) from 

— / (i) cos i] di (for«>l). 
Lc Jo 

Two-dimensional, undulated filaments can thus be mimicked 
by using Eq. (|3]l in which the amplitudes ajj are randomly cho- 
sen from a Gaussian distribution with mean value and stan- 
dard deviation (see Appendix B) 



The shape of cytoskeletal filaments such as F-actin is ther- 
mally perturbed by collisions with cytoplasmic molecules. 
The resistance of polymer chains to such thermal forces is 
characterized by their persistence length Lp, the characteris- 
tic arc length above which the filaments' tangent becomes un- 
corrected. Figure [T] shows a schematic of a snapshot of a 2D 
thermally undulated filament of contour length Lc, parame- 
terized by the arc length s and having a tangent angle 0{s). 



(4) 



Since dx(i) = cos[0(i)]di and dy{s) = sm[9{s)]ds, the fila- 
ment's coordinates x{s) and y{s) follow directly from 



x(s)^ / cos[0(/)]d/ 

Js'=0 



sin [e{s')]ds'. 



(5a) 



(5b) 




FIG. 1: Schematic of a thermally undulated filament of contour 
length Lq- The arc length along the filament is denoted by j e [0,Lc] 
and its tangent relative to the x-axis by the angle 6 (s) . The unit vec- 
tor along the end-to-end direction is mg. 



Finally, without loss of generality, we can require the end-to- 
end direction to coincide with the x-axis, i.e., y{Q) — y{Lc) = 
0. Application of the small-angle approximation in Eq. ( Bbl ). 
along with then yields 



9is')ds' = a/2L^«o =^ flo ~ 0. 



(6) 



Figure |2] shows a filament generated with the procedure de- 
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FIG. 2: A semiflexible filament, randomly generated using Eqs. 
JSjl with Lp = Lc = 10 /im, for different numbers of terms, A', in the 
Fourier series. 



scribed above using Lp = Lc ~ 10 /xm. The six curves all dis- 
play the filament in the same configurational state, i.e., with 
equal (random) amplitudes aj], but taking into account a vary- 
ing number of terms in Eq. As can be seen, the shape 
of the filament does not significantly change for A^ > 10. 

As a result of thermal undulations, the end-to-end distance 
ro will be smaller than the contour length Lq- To estimate 
the mean value of ro we can use the mean-squared end-to-end 
distance given by— 





■ s'-s' 


/ / exp 

Js=0 Js'=s 


2Lp 



- 4Lp ( 2 < exp 



2Lp 



- 1 



ds'ds 



Lc 
Lp 



(7) 



For filaments as in Fig. |2] having Lp = Lc ~ 10 jJ-m, (ro) 
is calculated using (|7]i to be about 9.23 fim. We have verified 
this by calculating ro = y/x^{s ^ Lc) + {s ^ Lc) for 10^ fil- 
aments generated as described above with A^ = 10, resulting in 
an avarage value of ro ~ 9.25 jJ.m^ 

The distance over which the ends of the filaments should 
be displaced to reach the straight configuration is Lc — ro, the 
slack distance ^ . The relative slack OJ, defined as 



(8) 



has an average value of 8.1 % for these filaments. To inves- 
tigate the dependence of dJ on the number of terms A^ in the 
Fourier series, we plot the relative slack probability density 
function p{{il) for A^ = 3, 10 and 100 in Fig. [3] It can be seen 
that the distribution shifts to larger CJ for larger values of A^. 
However, no significant difference is found between the dis- 
tribution functions for A^ =10 and 100. 



FIG. 3: Probability density function p{CS) of the relative slack C7 
in semiflexible filaments with Lp = Lq = 10 ^m, calculated using 
Eq. ([8) for 2 X 10"* filaments for A? = 3 (red), 10 (cyan) and 100 
(green). The black, solid curve shows the theoretical distribution 
function, according to Eqs.(|9) and ( llOt . 



One can also compare the relative slack distributions with 
the 2D radial distribution function given by Wilhelm and 
Frey^, 



^(ro;Lc;Lp) 



exp 



2Lf 



Lc ^ {2k 

k>0 



m 



1 



2kkl [2Lp(Lc-ro)]V4 



2Lp(Lc-ro) 



D 



3/2 



2{k+\)Lc 
yj2Lp{Lc- ro) 



(9) 



with Dt^i2{x) a parabolic cylinder function. A derivation 
of this equation is given in Appendix C. This function of 
ro = Lc I {I + (J3) , normalized as '^{ro;Lc',Lp)dro = 1, can 
be transformed into the distribution function for the relative 
slack: 



p{(D;Lc;Lp) 



Lc 



(1+0J)2 Vl+d 



Lc 



\Lc-M 



(10) 



The function /?(CI; 10 /im; 10 /im) is plotted in Fig. |3](solid, 
black curve) and agrees well with the generated distributions 
for large A^. 



III. THE STATIC AND DYNAMIC MECHANICAL 
BEHAVIOR OF SEMIFLEXIBLE FILAMENTS 

In this section we study the mechanical response of individual 
semiflexible filaments subjected to a tensile load. As men- 
tioned in Sect. I, this behavior depends on whether the un- 
dulation dynamics is taken into account. In addition, we in- 
vestigate the enthalpic contribution of stretching in extensible 
filaments. 
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We start by investigating the purely static^ mechanical be- 
havior of originally undulated, inextensible filaments. In the 
absence of an externally applied force, the chain in an undu- 
lated configuration is characterized by the set of amplitudes 
{flj]}, as descibed in Sect. II. Next, the chain is subjected to 
a tensile force /c along its end-to-end direction (the x-axis in 
Fig. lU, which induces a bending moment that tends to flat- 
ten out the filament. As a result, the configurational state will 
change to {an{fc)} with a decrease in the absolute value of 
each mode amplitude, i.e., |fl„(/c > 0)| < |flj]|. Describing the 
shape of the filament at each force level /c as a Fourier series 
according to Eq. (O, the transverse deflection y{s;fc) is given 
by 



2 y, a„(^fc) . 



(11) 



The dependence of each mode amplitude on the applied force 
can be calculated using the beam equation as follows. For an 
external force /c , the induced bending moment at a point s 
along the contour is given by fcy{s\fc), which has to be equal 
to the bending moment corresponding to the curvature of the 
filament, i.e.. 



fcy{s-Jc 



ds 



(12) 



Substitution of Eq. ^ (using a„{fc) instead of a^) and 
Eq. ( fTTT ) into Eq. ( fT2] i results in 



ifc + ml)an^O («>1). 



(13) 



By increasing the external force from fc to fc + d/c, mode 
amplitudes change to a„ +da„. We can then apply Eq. (flJl i 
to find a linear, first-order differential equation for each mode 
amplitude a„(/c): 



da„ 

a„ 



d/c 



(14) 



fc + ml ' 

With the initial condition fl„(/c = 0) = aj], the solution reads 



fln(/c) 



ml 



/c + ml 



(15) 



Equations ((Sj and ( fTSl l describe the shape of the filament at 
a given applied force, which enables us to calculate the force- 
dependent end-to-end distance. According to Eq. ( fSab . the 
end-to-end distance in the small-angle approximation is 



'"(/c) =x{s^Lc\fc)~Lc-- 



1 l-i-c 



which, after substitution of the Fourier series for 0, yields 

1 
2 



r{fc. 



(16) 



or, with the aid of (fTsT i, 



,0\2 



(/c 



mir^""^ 



(17) 
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FIG. 4: The force-extension relation of a randomly generated, in- 
extensible filament (dashed curve) and extensible filament (solid 
curve), with Lp = Lq = 10 and Kf = 8.53 x 10"'^ Nm^. The 
displacement is normalized with the filament's slack distance ^ 
(= 0.457 flm) and the force is normalized with a value of /c.max = 
0.01 N. The dotted line indicates the vertical asymptote for inex- 
tensible filaments. The value of the stretching stiffness used for the 
extensible filament is = 16 N. 



Alternatively, we can cast the force-extension relation m(/c) 
'"(/c) ^ '"0 = '"(/c) ^ '"(0) in the form 



E 

n>l 



f-+2f,ml ,n.2 
[U + mlY ^ 



(18) 



Figure |4] shows the mechanical response, fc{u), accord- 
ing to Eq. ( fTSl l of a randomly generated, inextensible fila- 
ment with Lp = Lc = 10 lim and K{ = 8.53 x 10"'^ Nm^ 
(dashed curve). The displacement u is normalized with the 
slack distance ^ = Lc — r{Q) = ^ L,i>i having a value 

of 0.457 jUm for this specific filament, and the force is nor- 
malized with /c.max = 0.01 N. Initially, undulations are being 
pulled out at relatively small forces, but as the displacement 
increases and the filament becomes (close to being) straight 
(i.e., M — > ^ or r ^ Lq), the filament locks and the force di- 
verges, since the inextensibility allows for no additional axial 
straining (see vertical, dotted Une). If, however, a chain has 
a finite stretching stiffness jU, it can be stretched beyond its 
slack distance ^ . In this case, the force-dependent end-to-end 
length r^(/c) can be constructed from Eq. (fTTI i for inextensi- 
ble filaments, by applying the transformation 



rpi{fc. 



1 



/c 



1 



(19) 



according to Storm and coworkers.— The resulting force- 
extension relation for /i = 16 N (see Appendix A) is plot- 
ted in Fig. |4] (solid line). For small displacements the me- 
chanical response of the extensible and inextensible filament 
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is identical, since this behavior is then dominated by internal 
bending. However, as the filament straightens, the force in- 
creases and the filament extends due to its nonzero compliance 
(i.e., finite stretching stiffness). For m/^ > 1.0, the mechan- 
ical response of the filament gradually becomes linear with a 
stiffness d/c/dw approaching ^/Lq^ As a check, the force- 
extension relation was also calculated numerically by dividing 
the filament into 500 Euler-Bernoulli beam elements of equal 
length (see reference— for details); the force-extension curves 
in Fig. ID were reproduced to within 2%. 

Now that we have determined the characteristics of the ini- 
tial shape of an individual filament and its mechanical re- 
sponse to an externally applied tensile force, we can deter- 
mine the ensemble-averaged force-extension relation of inex- 
tensible filaments. In the static approach adopted so far, the 
interaction of chains with the surrounding fluid is taken into 
account only via their initial, undulated configurational state, 
as can be seen from Eq. ( fT7] i. The chain's initial free energy 
is given by the Hamiltonian 



25 



20 



15 - 



10 - 



5 - 







1 


- 6000 




I 1 

1 1 
1 1 
1 1 


_ 4000 




1 1 
1 1 


- 9- 


' 1 


1 1 
; 1 


2000 


J 


/ / 
/ / 


- 



- 




/ / 


0.99 1 


/ / 
/ / 
/ / 
/ / 




Dynamic 


/ / 




/ 


/ Static 








, 1 


1 r r ~ "> ' ~'' ' 1 ' L. 






0.92 0.94 0.96 0.98 

(r)/L^ 



^„ 1 fdoV ^ 



(20) 



for the bending energy in the absence of tension (/c — 0). 
By substituting Eq. (|3]l into Eq. (l20l i and by setting each har- 
monic energy mode equal to k^T /2 (equipartition theorem), 
the mean-squared value of becomes 



FIG. 5: The ensemble-averaged mechanical response of inextensi- 
ble semiflexible filaments having Lp = Lc = 10/im and Kf = 8.53 x 
10^'^ Nm^ in terms of the end-to-end distance (r) versus the dimen- 
sionless force (p = fcL^/ {KfTt~). Results are plotted for the static 
model, using Eq. ([23} [solid curve], and for the dynamic model us- 
ing Eq. J30t [dashed curve]. The inset shows the corresponding aver- 
age mechanical response of extensible semiflexible filaments, using 
Eq. (T9I1 with /X = 16 N. 



1 



(21) 



This result can be substituted into Eq. ( [TtI i. yielding the force- 
dependent, ensemble-averaged end-to-end distance 



{r) i(p)=Lc- 



n>\ 



(22) 

where (p is the force normalized with the critical Euler buck- 
ling force, i.e., (p = fcL^/ (KfTi^). The summation in Eq 
up to ^ 00 can be carried out analytically, giving 



Second, the initial value of the axial stiffness Gb = d/c/d(r) 
resulting from internal bending by pulling out the thermal un- 
dulations is found to be 



G« = Gb((r) = (ro)) = 



d(p 



d(p d/c 



90LpKf 



(25) 



/c=0 



The diverging behavior for large forces can be found from 
Eq. ( l23b . Since coth[7ry/^] 1 and sinh"^[7ry/^] for 
large (p, the limiting behavior can be expressed as 



, , , X /coth[7r,/ffll „\ 

Figure|5]shows the static, mechanical behavior according to 
Eq. ( l23T l [solid curve]. The response is linear at small forces 
and diverges as the average end-to-end distance approaches 
the contour length. The linear response at small forces can be 
found by a Taylor expansion of Eq. ( [23] l, 



Lc 



12Lf 



L^n (p 
90Lp ' 



(24) 



First, we conclude from this, or directly from (r) (0) in 
Eq. ( l22b . that the average slack distance for 2D filaments is 
((^) = L^/ (12Lp). This value is also obtained by expanding 
the square root of the result in Eq. (|7]i in a Taylor series.— For 
Lc ~ Lp ~ 10 /xm the initial, average end-to-end distance is 
9.17 jxm, close to the value of 9.25 jim we found in Sect II. 



(r)((p)^Lc- 



(26) 



implying that the force-extension relation diverges as fc oc 
{Lc — {f))^^ as (r) approaches Lc, just like the worm-like 
chain model.-- 

These results can be directly compared with the results 
derived by MacKintosh and coworkersji^ In their entropic 
model, filaments continually undergo thermal bending mo- 
tions as they are subjected to an external tensile force. At 
each instance and applied force the configurational state of 
the filament changes, as a result of the interaction with the 
surrounding fluid. In this dynamic model, each mode ampli- 
tude changes in time and with force, i.e., a,, = an(f;/c). The 
free-energy functional of the filament is written asi^ 



e^{s)ds, (27) 



6 



where the first term is the filament's bending energy and the 
second term is the work of contracting against the appUed ten- 
sion. Using this Hamihonian, the force-dependence of the 
mean-squared value of a„ can be found from equipartition at 
each force, and reads 



k')(/c) 



1 



1 



Lpql+fc/Kf' 



(28) 



Then, by making use of Eq. ( fT6b . the time- or ensemble- 
averaged end-to-end distance (r) becomes^ 



(29) 



and can be rewritten in the following form: 



{f){(p)=Lc- 



Ll 



4-Lpn^(p 



(7rV?coth[7rV^ - 1). (30) 



The dynamic behavior, shown in Fig.|5]by the dashed curve, 
is thus quite similar, but not equal to the behavior according 
to the static description. The small-force response can be ap- 
proximated by 



(<?) 



12Lp 



180Lp ' 



(31) 



so that the slack distance in the dynamic model is equal to that 
in the static model, cf. Eq. ( l24b . However, the initial stiffness 



ISOLpK-f 



(32) 



is a factor two larger than the static stiffness, cf. Eq. dZSl ). 
At large forces, the diverging response is also similar to what 
we found in the static model, but from Eq. ( l30l l the limiting 
behavior near full extension. 



(r)((p)^Lc- 



(33) 



has a stiffness that is a factor four higher than in ( |26] |. The 
dynamic-to-static stiffness ratio, G/G, is shown as a function 
of (r) /Lq in Fig. |6] (solid curve). As mentioned above, the 
ratio gradually increases from 2 at (r) /Lc = 1 — ic/(12Lp) 
to 4 in the Hmit (r) /Lq 1. 

The equality of the slack distance in both models is sim- 
ply due to the fact that, in the absence of an external force, 
the Hamiltonian used in the dynamic model, Eq. dZTT i. reduces 
to the bending energy, Eq. (|20] |. that is also used for equipar- 
tition in the static model. However, as soon as the chain is 
subjected to a tensile force (no matter how small in magni- 
tude), the static and dynamic mechanical responses deviate. 
This deviation is caused by the difference in the chain's en- 
ergy used for equipartition. In the static model, only the initial 
undulations are taken into account by equipartition at /c = 0. 
Next, each frozen-in filament is subjected to a tensile force, 
thereby reducing the internal bending energy of the system. 
In the dynamic model, equipartition is imposed on the total 
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FIG. 6: The the evolution of the dynamic-to-static stiffness ratio 
(G/G) and external work ratio (ly/W) with(r) /Lq- Results are 
plotted for inextensible filaments (G/G solid curve; W /W dashed 
curve) and extensible filaments (G/G dotted curve; W /W dash- 
dotted curve). The parameters used for the calculation: Lp = Lc = 
10 lim, Ki = 8.53 X 10"'^ Nm^ and /i = 16 N. 



free-energy functional, comprising both the bending energy 
and the chain's work of contracting against the applied force: 
on average each harmonic mode in the chain's total internal 
energy remains kgT /2 during pulling, even though its bend- 
ing component reduces. At each force, there must therefore be 
sufficient time for a complete energy redistribution between 
different bending modes and between the bending and tension 
energies, to ensure that the chain's energy [Eq. dZTl iI remains 
constant. Compared to static filaments, a larger force is re- 
quired to reach a certain end-to-end distance in order to over- 
come the effect of dynamic undulations, thus resulting in a 
higher stiffness. 

Since energy is the origin of the difference in static and 
dynamic mechanical response, it is of interest to study the ex- 
ternal work W needed to reach an average end-to-end distance 
(r), which is calculated as 



W{{r)) 



r) 
ro) 



(p{s)ds 



(34) 



with (ro) = Lc — Lq/ (\2Lp). By numerically inverting 
Eqs. ( |23] |, ( |30] | and by integration with Eq. (|34] |. it can be 
shown that W{{r)) diverges as (r) approaches Lq, both in 
the static and in the dynamic model. Therefore, the dynamic- 
to-static work ratio, W /W, approaches the same value as the 
stiffness ratio, i.e., 

lim {W/W)= lim (G/G) =4. 

The dependence of the work ratio on (r) for inextensible fila- 
ments is shown in Fig. |6] (dashed curve) and indeed increases 
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from 2 at (r) = (ro) to 4 near full extension, similar to the 
stiffness ratio. Compared to static chains, the energy needed 
to extend dynamically undulating chains close to their contour 
length is about four times larger. In both cases, however, the 
completely straight configuration will never be reached as a 
result of the energy divergence. 

The origin of this divergence is different in the two models 
and can be explained by allowing thermal bending motions to 
excite only the first order mode ai in the Fourier series Q. 
The average value of the end-to-end distance in the absence 
of force then is (ro) = Lc — L^/{2k^Lp), which follows di- 
rectly from Eqs. ( |22] | and ( |29] l. The force-extension relations 
can be inverted analytically and the external work according 
to Eq. ( l34b can be readily calculated. For static chains, the 
energy needed for complete stretching is 




and, thus, is equal to the average initial energy stored in bend- 
ing. For a finite number of allowed modes, static chains can 
be stretched to their contour length by an energy input equal 
to the initially stored bending energy. The origin of the en- 
ergy divergence for static chains is merely the result of allow- 
ing an infinitely large number of bending configurations in the 
Fourier series (N °o). However, for dynamic chains the en- 
ergy needed for full stretching is 



lim 



-1 



In^LpiLc-s) 



di, 



and diverges. This is a result of the fact that equipartition is 
applied on the total free-energy functional [Eq. (|27]|1. Hence, 
dynamic chains can never be pulled straight completely, even 
if only a single bending mode is allowed. 

In case of extensible filaments we have to incorporate the 
axial stretching energy 



2 70 



fALr 



ds, 



(35) 



where d/ / As is the relative length change along the filament 
and ALc the total increase in contour length. Equation (l35T l de- 
termines the Hookean response of the filament with a stretch- 
ing stiffness [i. The ensemble-averaged mechanical behav- 
ior of such filaments can be found by applying Eq. ( fT9] l to 
Eq. (|22] | for static chains and to Eq. ( |29] l for dynamic chains. 
The results are shown in the inset of Fig.|5] with = 16 N. As 
(r) approaches Lc, the mechanical response is clearly domi- 
nated by axial stretching of the filaments. The energy stored 
in stretching dominates the total energy and the deviation be- 
tween the dynamic and static model vanishes around (r) = Lc- 
This can be shown even better by studying the dynamic-to- 
static stiffness and work ratios for extensible filaments, that 
are also shown in Fig. |6] (dotted and dash-dotted curves, re- 
spectively). While following the curves for inextensible fila- 
ments at small end-to-end displacements (bending-dominated 
regime), both ratios decrease to one around (r) = Lc, where 
the stretching-dominated regime starts. 



The stiffening of individual filaments cannot be easily over- 
estimated: the average axial stiffness increases from Gb = 
7.7 N/m for static filaments, Eq. dHll, or Gb = 15.4 N/m 
for dynamic filaments, Eq. ( l32b . in the bending-dominated 
regime, to a value of Gs = [ijLc = 1-6 x 10^ N/m in the 
stretching-dominated regime. Since Gb is six orders of mag- 
nitude smaller than Gs, we will neglect from now on the en- 
tropic, axial stiffness due to internal bending. For an individ- 
ual filament with a slack distance ^, the dependence of the 
axial stiffness on the end-to-end displacement u is thus ap- 
proximated by 



for u <^ 
II /Lc for M > 



(36) 



Since the probability of a filament having a slack distance 
between ^ and ^ +d^ is equal to ^{Lq — ^;Lc;Lp)d^ [cf 
Eq. (|9]l], the average stiffness of an ensemble of chains with 
different slacks can be expressed by 

(Gi)(«;Lc;Lp)= r'^(Lc-^;Lc;Lp)Gi(M,<^;Lc)d^ 
Jo 



^ /V(Lc-^;Lc;Lp)d^ (37) 
Lc Jo 



Lc 

This result can be integrated to yield the average force- 
extension relation: 

(/,)(m;Lc;Lp)= r (Gi) (M';Lc;-Lp)dM' (38) 

Ju'=0 

It is important to note the difference in averaging procedures 
used above: in Eqs. ( |22] | and ( |29] l the average is taken over 
all end-to-end lengths at constant force, whereas in Eqs. dJTl i. 
(|38T l the average is taken over all forces at constant displace- 
ment. 

The average stiffness, (Gi) (m; 10 jim; 10 jJ-m) according to 
Eq. (|37] | is plotted in Fig.|2ta) [solid circles]. As can be seen, 
the average stiffness gradually increases from for small dis- 
placements to a saturation value of /i/Lc for large displace- 
ments. This enthalpic stiffening thus results from an increase 
in the number of filaments that are stretched beyond their 
slack distance. The corresponding average force-extension re- 
lation is shown in Fig.|7tb) [solid circles]. 

To investigate and check the mechanical response found in 
the analytical model, we have performed 2D finite-element 
(EE) calculations on individual filaments (/x = 16 N, iff = 
8.53 X 10-'^ Nm^). Filaments were generated as described 
in Section II using Eqs. (O and ^ with Lp = Lc = 10 jUm, 
fixed at one end and pulled by prescribing the displacement 
of the other end along the direction of the initial end-to-end 
vector. The force-extension relation was calculated for 2000 
randomly generated filaments, from which the ensemble- 
averaged response was determined by averaging the force 
at constant displacement; the stiffness was computed from 
d (/( ) /du. The results are shown by the solid curves in Fig.|7] 
As can be seen, the FE-predictions correspond very well with 
the analytical model. What cannot be seen, because of the 
plotting scale, is that the initial average stiffness is 10^ N/m 
instead of the value 7.7 N/m calculated from Eq. ( |25] ). This 
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FIG. 7: (a) The ensemble-averaged stiffness {G\) {u;Lc',Lp) of in- 
dividual semiflexible filaments with Lp = Lc = 10 /im, and (b), 
the corresponding average force-extension relation (fc) {u;Lc\Lp). 
Results are shown of the analytical calculations using Eqs. ( I37t . 
J38b (closed circles), and from the finite-element (FE) calculations 
(drawn curve). The mechanical parameters used are: /i = 16 N and 
K{ = 8.53 X 10^'^ Nm^ (the latter value is needed in the FE calcula- 
tion only). 



is due to the difference in averaging procedure as explained 
above; if we average over displacement at constant force, the 
result for small displacements indeed corresponds to Eq. dZSt . 



IV. AFFINE MODEL FOR ELASTICITY IN 
SEMIFLEXIBLE BIOPOLYMER NETWORKS 

The ensemble-averaged force-extension relation of individ- 
ual filaments derived in Sect. Ill, enables us to calculate the 
mechanical response of cross-linked networks of such fila- 
ments under simple shear To do so, we follow the procedure 
used by Wu and Van der Giessen for the elasticity of rubber 
networks 

We consider a square unit cell of dimension W spanned by 
base vectors ei and 62 as indicated in Fig. [Sa). In the ini- 



arctan(r) 




FIG. 8: (a) Schematic representation of affine deformation of a two- 
dimensional unit cell under simple shear at a strain F. The initial 
end-to-end vector ro at an angle <I> with ej has transformed to a vector 
r = Fro at an angle (j). (b) An (undulated) filament with end-to-end 
length r at an angle (p under a tensile force fc in the deformed state. 
The area rw covered by the filament is, on average, constant under 
the assumption that the network is incompressible. 



tial state, the unit cell contains A^f filaments of length Lq ~ 
Lp = 10 /im, at random orientations. The areal and line den- 
sity of filaments in the network are defined as n ^Nf/W^ and 
p = nLc, respectively. The angle between the initial filament's 
end-to end-vector ro = romo and ei is denoted by <I>. 

Next, the cell is subjected to a simple shear of strain F, 
by fixing the lower boundary and displacing the top boundary 
over a distance WT, see Fig. [HJa). This deformation process 
is represented by the deformation gradient tensor F with com- 
ponents 



[Fij] 



1 F 
1 



(39) 



Note that the shear deformation process satisfies detF = 1 , ex- 
pressing incompressibility of the filamental network42 When 
it is assumed that the filaments deform in an affine manner, 
each filament's end-to-end vector Tq transforms to the vec- 
tor r = rm according to r = Fro, as sketched in Fig. [8ja). 
Filaments thus rotate from <I> to and undergo a stretch of 
A = r/ro that can be obtained from F through 



^m(FF'^ 



m, 



(40) 



where m = r/r = cos 0ei + sin ^62 is the unit vector along the 
filament's end-to-end direction in the current, deformed state. 
From Eqs. ( [39] l and ( l40b the stretch X can be expresssed in 
terms of and F as 



1 + tan2 (j) 



2Ftan^ + (l+F2)tan2 0' 



(41) 
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FIG. 10: Polar plot of the chain orientation distribution function 
FIG. 9: The stretch X = r/ro of a filament as a function of its orien- C((|);r) [CODF], at r = 0.0 (solid curve) and T = 0.3 (dashed curve), 
tation (j) e [-7r/2, 7t/2] in the deformed state at various levels of the The dotted lines indicate the principal directions at T = 0.3. 
shear T according to Eq. l l41t for affine deformation. 



A filament that is at an angle <j) in the current state of shear 
strain F is thus subjected to a stretch A that is given by 
Eq. Figure|9]shows the dependence of A on for various 
strains F. 

During shear, filaments also rotate and thus change their 
orientation. If C(0;F)d0 is the probability of a filament with 
an orientation between (j) and (j) +d(j) in the deformed state at 
strain F, it can be shown that C(0;F) follows directly from the 
stretch A by^ 

C((/);F)=CoA2(,/);F), (42) 

where the normaUzation value Co = l/nis the initial, uniform 
probability distribution. This probability function C(0;F) is 
also known as the "chain orientation distribution function" 
(CODF).i22i^ Figure [To] shows a polar plot of the CODF for 
F = and 0.30 in the angular range [Q,2n] .— For F = there 
is no deformation and filaments are randomly distributed in 
the network, indicated by the spherical distribution. As the 
network deforms, filaments rotate in the straining direction, 
resulting in a nonspherical distribution, its anisotropy grow- 
ing with increasing strain. For F = 0.3 the principal directions 
for which a minumum and maximum stretch is obtained, are 
indicated by the mutually orthogonal, straight dashed lines. 
For this strain, filaments at an angle of (p = 40.7°undergo a 
maximum stretch of A = 1.16 and filaments at an angle of 
(p = —49.3° a maximum compression (A — 0.86). 

Under the assumption of affine deformation, the Eqs. 
( 136137138141142b can be used to calculate the network re- 
sponse to simple shear We start by considering an individ- 
ual filament at an angle (p subjected to a stretch A (^,F) in the 
deformed unit cell of strain F, as sketched in Fig. ^h). As- 
suming the network is incompressible, the area per chain, rw, 
remains constant, and is thus equal to l/n. The force acting 



on the filament, /c, can be tranformed to the Cauchy traction 
(7c acting on the continuum in which the filament is embed- 
ded, C7c = fc/{bw), where w is the cross-sectional width and 
b is the unit out-of-plane length the force acts on. Similar 
to the approach by Wu and Van der Giessen^ we define the 
micro-stress tensor C7c by 

CTc = o-c(m®m)-pI, (43) 

with Gc = nA(0;F)/cro^"' since r ~ Aro. In Eq. ( |43] | the hy- 
drostatic stress p is included to account for the incompressibil- 
ity of the network. The micro-stress tensor is the contribution 
of a single filament oriented along in(0) to the stress of the 
network. Finally, with the areal density of chains dn having 
an orientation between ^ and + d^ given by 

d« = «C(0;F)d0, (44) 

the overall or macro-stress tensor O of the network can be 
evaluated by the average of the micro-stress tensor over 
the individual chains in the network, i.e., 

CT = - f Ocdn- (45) 
n J 

It is important to note here that since the force acting 
on an individual filament depends on its end-to-end dis- 
placement, it therefore depends on its current orientation (p 
and applied shear strain F through /c(0;F) = /c(m(0;F)) ~ 
/c(ro(A((/);F)-l)). 

We assume that the filaments in the network only 
stretch and compress upon mechanical loading, but do not 
bend/buckle. In Sect. Ill we showed that stretching by in- 
ternal bending can also be neglected, and that axial stretching 
with stiffness fi/Lc dominates the ensemble-averaged single- 
filament force-extension relation. The angular region in which 
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FIG. 11: (a) Shear stress T versus shear strain F response for biopoly- 
mer networks comprising filaments with Lp = Lq = 10 fim, Kf = 
8.53 X 10^'^ Nm^ and = 16 N, and (b), the corresponding shear 
stiffness G = dt/dF. Results are plotted for three filament densities: 
p=13 (black curves), 25 (blue) and 38 (red). The solid curves display 
the results from the discrete-network (DN) calculations, the dashed 
curves are the results from the affine-network (AN) model, Eq. i46) . 



Stretching occurs (A > 1) follows directly from Eq. dTTT i and 
is given by S [0,arctan(2r^')]. In this region filaments 
experience an elongation of their end-to-end distance; fila- 
ments outside this region are under compression {X < 1) and 
deform by soft (internal) bending and contribute negligibly 
to the overall network stiffness. Hence, we adopt the force- 
extension relation based on Eqs. ([36]l-(l38]l in the angular re- 
gion £ [0,arctan(2r^^)]. The dimensionless shear stress, 
T = GnbLc/jj., from Eqs. (|42]|-(|45]| can then be written as 

arctan(2r^') 

nil Lc J 



(46) 

where p = pLc = Ni{Lc/W)'^ is a dimensionless filament 
density. Figure fTTT a) presents the calculated T versus F-curves 
for three densities p :=:13, 25 and 38 (dashed curves). The cor- 



responding shear stiffness, G = dr/dF, is shown in Fig. fTTT b). 
As can be seen, stiffening of the network is indeed observed, 
resulting from the average stiffening of individual filaments 
in the network (Fig. |7]i, or said differently, stiffening results 
from an increase in the fraction of filaments that are stretched 
beyond their slack distance. 

These results can be directly compared to 2D numerical 
discrete-network (DN) calculations using the finite-element 
method.^ Filaments of contour length Lc and persistence 
length Lp are randomly placed into a square cell of dimen- 
sion W at random orientations, with proper account of period- 
icity. Each filament is discretized into 15 equal-sized, Euler- 
BernouIIi beam elements that account for stretching (with ax- 
ial stiffness jj.) and bending (with stiffness fCf). In accordance 
with Sect. II, 10 Fourier modes are used to describe the ini- 
tial shape of the undulated filaments. The points were fil- 
aments overlap are considered to be stiff cross-links, where 
the rotation and the displacement of the two filaments at the 
cross-link is equal. Figure [T2l' a) shows an example of a ran- 
domly generated network of density p = 1 3 that is in its initial, 
stress-free configuration. Filaments crossing either the upper 
or lower boundary of the cell are perfectly bonded to rigid 
top and bottom plates, respectively. Next, the top plate is dis- 
placed horizontally relative to the bottom plate over a distance 
rw, corresponding to an applied shear strain of F. Geometry 
changes are accounted for by an updated Lagrangian finite- 
strain formulation. Figures [TSlb) and (c) show the network 
of Fig. [TSla) in the deformed state at shears of r=0.125 and 
0.285, respectively. The macroscopic shear stress T is calcu- 
lated from the cumulative, horizontal reaction force of nodes 
at the top plate, divided by the cell width W. Convergence 
studies were performed to ensure that the results are not af- 
fected by the number of elements per filament and the cell 
size W. 

The stress-strain response (averaged over 10 different, 
random realizations) at three network filament densities of 
p =13, 25 and 38 is included in Figure [TT] (solid curves). It 
can be seen that, at a relatively small density (p =13), the 
discrete-network (DN) calculations result in a much softer 
response compared to the affine-network (AN) calculations. 
Closer examination of the DN results reveals, first, that the 
small-strain mechanical response is dominated by bending 
and buckling of filaments^^ thus yielding a soft response. 
Second, at somewhat larger strains, filaments perform addi- 
tional, nonaffine motions during shearing: filaments reori- 
ent themselves in the direction of straining by rotations and 
translations. These nonaffine, geometrical network reorien- 
tations can be clearly seen by comparison of Fig. [TSlb) with 
Fig. [12 a). Both effects, bending/buckling of filaments and 
local network reorientations, are not taken into account in the 
AN model, where instead, at small strains, all filaments with 
an orientation between and 7r/2 are subjected to stretch. 
The response in the AN model is therefore dominated by 
stretching of filaments, contrary to what is observed in the 
DN computations. However, as the strain increases, more fil- 
aments become oriented in the direction of straining, as seen 
in Fig. [T2l c): strings of stretched filaments connect the top 
and bottom plate of the cell. Thus, at larger strains, the re- 
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FIG. 12: Shear deformation of the discrete network (DN). (a) Initial, F = 0; (b) intermediate, F = 0. 125 and (c) large strain, F = 0.285, network 
configurations for a randomly generated network of density p = 13, with W = 4Lc = 40 flm. 



sponse according to the DN calculations is also dominated 
by stretching of filaments, which explains the stiffening in 
the DN curves in Fig.[TT|at low densities. Since the relative 
amount of filaments in stretching remains small compared to 
that of the AN model, the stiffness determined from the DN 
calculations remains small compared to the stiffness from the 
AN model (black curves in Fig.fTTI). 

In contrast to small densities, at larger densities (e.g., p = 
38) the DN calculations exhibit a higher stiffness compared to 
the AN calculations (red curves in Fig.fTTTi. The reason for this 
is that the cross-link density increases with filament density, 
with two effects. First, the higher cross-link density results in 




FIG. 13: The product pic versus the filament network density p (in 
/xm^'), determined from and averaged over randomly generated net- 
works (solid squares). The dashed line indicates the limiting, high- 
density value of 7t/2. The inset shows a random network at a density 
of p = 20 /im" 1 (p = 200), with W=Lc = WfXm. 



a more rigid network structure in which each filament is cross- 
linked to more filaments, thereby decreasing the freedom of 
filaments to undergo nonaffine reorientations. The transition 
from a bending-dominated regime to a stretching-dominated 
regime therefore shifts to smaller strains as can be seen by 
comparing the solid DN curves in Fig.[TT] Second, chain seg- 
ments between cross-links contain less slack than the original 
filaments: segments are thus stretched beyond their slack dis- 
tance at smaller strains leading to a higher stiffness. In the AN 
model, this reduction in slack is not taken into account, but 
cross-links are assumed to be at the ends of the filament, in- 
dependent of the density of filaments in the network. Density 
enters the AN model only through the prefactor in Eq. ( l46b . 
which explains the (self-) similar shape of the AN curves in 
Fig.E] 

To correct for the decrease in slack in the AN calcula- 
tions, we first have to determine the average cross-link dis- 
tance, Ic, in the network as a function of the density p. 
The values obtained from a series of generated networks are 
shown in Fig. [13] as the product pk versus p. The product 
gradually increases from 1.2950 ±0.0325 at p = 1.25 /im^' 
(p = 12.5) to 1.5618 ±0.0034 at p = 35.0 ixm-^ (p = 350) 
and thus approaches the limiting, 2D theoretical value of 
^^2.iMi The inset shows an undeformed network with W = 
Lc = 10 jj.m at a high density of p = 20.0 /xm^', exem- 
plifying the short average distance between the cross-links 
(l^ = 0.0775 ± 0.0003 ^m < Lc). Next, the distribution in 
end-to-end lenghts ro of segments at a density p is given 
by ?^(ro;/c(p);ip) from Eq. with Lc replaced by kip)- 
From this, the ensemble-averaged mechanical behavior of a 
segment of mean length lc{p), i.e., (d) (m;/c(p);Lp) and 
ifc) (M;/c(p);ip), follows directly from Eqs. This 
then serves as input for the density-corrected AN response 1 
given by Eq. (|46]|.— The corrected AN response is plotted 
in Fig. [T4]for three densities p = 12.5, 25 and 37.5 (dashed 
curves). For comparison, the figure also displays the DN cal- 
culations from Fig.[TT](solid curves). As can be directly seen, 
all AN curves have gone up due the reduction in slack: already 
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FIG. 15: Initial, small-strain stiffness Go versus the density p, cal- 
culated using the DN model (solid squares). The values of Go result 
from an average over 10 random networks in case p < 100 and 5 net- 
works for higher densities. The dashed line is the affine upper limit 
according to Eq. J48t 



as in tension, i.e., /i/Zc, corresponding to an axial force of 
(/i//c)M. Hence, the integration should be performed over the 
entire angular range, [0, tt], and with the segmental force equal 
to /x[A(0;r) — 1], the network response is evaluated as 



FIG. 14: (a) Shear stress T versus shear strain F response for biopoly- 
mer networks comprising filaments with Lp = Lc = 10 ]lra, Kf = 
8.53 X 10^^^ Nm^ and ^ = 16 N, and (b), the corresponding shear 
stiffness G = dr/dP. Results are plotted for three filament densities: 
p=12.5 (black curves), 25 (blue) and 37.5 (red). The solid curves 
display the results from the DN calculations (same as in Fig.ll It. the 
dashed curves are the results from the density-corrected AN model. 



at small strains filamental segments are stretched beyond their 
slack distance resulting in a higher stress and stiffness. Af- 
ter the density correction, all AN calculations exhibit a higher 
stiffness than the DN calculations, irrespective of the density: 
as said before, this can be attributed to nonaffine network de- 
formations (bending/buckling, local network reorientations) 
in the DN calculations that are not taken into account in the 
AN model. 

In calculating the density-corrected response, the angular 
integration was performed in the interval [0,arctan(2r^')], 
see Eq. (l46T l, assuming that segments still contain a consid- 
erable amount of slack, resulting in a very compliant response 
in case such a segment is under axial compression. How- 
ever, as discussed above, when the filament density increases, 
segments between cross-links become shorter and virtually 
straight. For a network consisting of only straight filaments, 
segments under compression exhibit the same axial stiffness 



n 

T=- /A-^(0;r)(A((/);r)-l)sin(/)cos^d0. (47) 

7t J 



Eq. ( |47] | thus gives the limiting, affine mechanical behav- 
ior of biopolymer networks of high density. Whether or not 
filaments are undulated and contain slack is not relevant at 
this point, since the segments of average length lc ^ Lc are 
straight at high densities [see inset Fig. (fT3bl. As we showed 
before, the DN calculations exhibited nonaffine deformation 
characteristics at low densities. However, at higher densities 
the DN calculations are also expected to show an affine be- 
havior To check this, we have investigated the initial, small- 
strain network stiffness. Go, as a function the density p, see 
Fig. [15] (solid squares). The affine value follows directly from 
the derivative of Eq. (l47b by a first-order Taylor expansion of 
the integrand with respect to F, and equals^ 

Go,aff=|. (48) 

Note that, using the high-density limit pk — n/2r^the stiff- 
ness can be rewritten as Go.aff = ?fjti/(16/c), corresponding to 
the result found by Head et al. The normalized affine limit, 
Go.aff. is included in Fig.[T5]as the dashed line. It can be seen 
that, as the density increases, the DN calculations indeed ap- 
proach the AN limit. At a density of p = 350 the DN calcula- 



13 



tions result in a stiffness of Go = 38 .0 ± 0.6, close to the affine 
value of 43 .75. 

Finally, it is interesting to see that in the networks' stretch- 
ing regime, i.e., at large strains, the network stiffness (both 
AN and DN) continues to increase despite the fact that the 
individual segments have a constant axial stiffness jx/lc in 
stretching [see Fig. [T4T b)1. From this one would expect the 
stiffness to converge to a fixed 'steady-state' value. However, 
due to nonlinear geometrical effects at large strains, the stiff- 
ness increases with straining according to 

Gocp r2(i+r2)-V2^ (49) 

as derived in Appendix D. From further investigation of 
Eq. (l47T i it follows that the stiffness starts to level off at very 
large, experimentally irrelevant strains of F « 1.5. 



V. CONCLUDING REMARKS 

We have performed a detailed investigation of the mechan- 
ical stiffening behavior of 2D semiflexible biopolymer fila- 
ments and networks of such filaments under simple shear 
Such polymers undergo thermally excited bending motions, 
brought about by collisions with molecules of the surround- 
ing fluid in which the filament is immersed. For inextensible 
filaments subjected to a tensile force, these undulations can 
be pulled out at the cost of an external energy input, result- 
ing in an axial stiffness that can be calculated for two specific 
cases. First, we can assume no undulation dynamics to take 
place during tension. In this case, the filament is isolated from 
its surrounding fluid and the random, undulated state present 
in the filament is subsequently pulled out when subjected to a 
tensile force. In this static case, the axial stiffness results from 
internal bending with bending stiffness Kf; using equilibrium 
mechanics, the force-extension relation can be calculated and 
averaged over an ensemble of filaments. Second, undulation 
dynamics can be taken into account during pulling. By apply- 
ing equipartition on the filament's total energy functional at 
each force, the ensemble-averaged, dynamic force-extension 
relation can be derived. This so-called entropic stiffening was 
previously derived by MacKintosh and coworkers 

The resulting average mechanical response of single fila- 
ments in these two cases turns out to be very similar At 
small forces and end-to-end displacements, both responses 
scale with K^, with the initial dynamic stiffness being a factor 
2 larger than the static one. Close to full stretching, as the av- 
erage end-to-end length (r) approaches the contour length Lq, 
the force diverges as /c <^ (Lc ~ ('") ) ^ in both models, charac- 
teristic for the worm-like chain model, with the dynamic/static 
stiffness ratio approaching 4. This indicates that the static and 
dynamic approaches are qualitatively very similar, capturing 
the same physical dependence on the mechanical and geomet- 
rical parameters Kf, Lp and Lq, but differ quantitatively by a 
factor 2 to 4. 

In case of extensible filaments of stretching stiffness ji, en- 
thalpic axial elongation dominates the filaments' mechanical 
behavior close to full stretching and both models thus exhibit a 



quantitatively identical mechanical response. In addition, we 
have shown that for 2D filaments, one can very well describe 
the ensemble-averaged force-extension relation by completely 
neglecting contributions from internal bending and by merely 
taking into account the axial stretching of filaments and the 
distribution in slack. 

Two-dimensional, cross-linked networks that have these 
semiflexible filaments as building blocks stiffen under an ap- 
plied shear strain. Often, stiffening is attributed to the en- 
tropic stiffening of individual filaments undergoing affine de- 
formations. For the 2D networks under investigation, we con- 
clude by comparing analytical, affine network (AN) calcula- 
tions with discrete (finite-element) network (DN) calculations 
that this cannot be the case. For networks of low filament 
density, the mean cross-link distance lc is relatively large and 
the resulting segments between cross-links contain a relatively 
large amount of slack that can be pulled out. In this case, one 
could argue that entropic stiffening would be the main ori- 
gin of stiffening. However, the DN calculations show that at 
these densities stiffening lies in the network rather than in its 
constituents: at small strains the mechanical response is dom- 
inated by bending/buckling of filaments and nonaffine rota- 
tions and translations. These local network rearrangements 
induce a transition from a bending-dominated regime at small 
strains to a stretching-dominated regime at larger strains in 
which strings of connected filaments dominate the mechanical 
response. At higher network densities, the number of cross- 
links increases and the mean cross-link distance lc decreases, 
so that the shorter segments between cross-links contain less 
slack. The higher network density then results in a mechani- 
cal response that is more affine in nature, with the mentioned 
reduction in slack increasing the enthalpic stretching contri- 
bution of the segments. By accounting for the reduced slack 
in the AN calculations, the DN results were severely over- 
estimated, showing that, for the densities investigated, the 
discrete-network architecture plays a key role. Only in the 
limit of extremely high densities do the DN and reduced-slack 
AN calculations coincide. 

Hence, we conclude that, for physiological relevant densi- 
ties, stiffening in 2D networks results from nonlinearities in 
the network response, rather than from entropic stiffening of 
individual filaments. Note that these conclusions might be 
different for three-dimensional networks, since the amount of 
slack in filaments is significantly higher than in 2D. Current 
work focuses on this issue in 3D networks. 
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Appendix A 



Appendix C 



While real semiflexible filaments are, evidently, three- 
dimensional objects that are often idealized as rods of diame- 
ter d, 2D models require some 'projected' thickness t for the 
specification of the axial and the flexural stiffness in terms of 
Young's modulus E, the area of a filament's cross section. A, 
and its second moment of inertia, /. This 'projection', how- 
ever, is not trivial nor unique. While the moment of inertia of 
a circular cross-section of diameter d is I = 7rc/'*/64, that of 
a projected filament with thickness t is I = bt^ / 12 with some 
out-of-plane width b. Aiming at actin filaments which typi- 
cally have a diameter of c/ = 8 nm, we use 2D filaments with 
f = 8 nm. When a typical value of £ = 2 GPa is used for the 
Young's modulus of actin, the bending stiffness Kf of such 2D 
filaments with a unit out-of-plane thickness (i.e., = 1 m) is 



bP 



8.53 X 10"'^ Nm- 



(50) 



If we would use this result in Eq. (|2|i, the persistence length 
would be 2.1 x 10'* m, much larger than the contour length 
Lq and filaments would simply be straight, which they are not 
in 3D. Therefore, the persistence length is scaled back to the 
contour length (scale factor: 4.8 x 10^'") to generate semi- 
flexible filaments for which Lp = Lq = 10 jim, as described 
in section 11. 

For the calculations performed in sections 11, 111, and IV 
we use the same input value of Lp = 10 /xm to generate the 
undulated filaments (i.e., the initial configurational state {a^,}, 
determined by this value of Lp) , but for the bending stiffness 
Kf in the beam-equation [Eq. ( fT2b l we use the value calculated 
in Eq. ( |50l ). For this reason Lp and Kf are treated as separate 
parameters in the eqs. (l2ni-(|26]l. With the same 2D repre- 
sentation as underlying Eq. dSOl l. the value of the stretching 
stiffness EA is 



li=EA = EbtKi 16 N. 



Appendix B 



(51) 



For a Gaussian distribution with mean (aJi) ~ and standard 
deviation s„ , the probability p{a^)da'^ of finding a mode am- 
plitude between aj] and aj] + da^ is 



p(a«)d«° 



: exp 



The mean-squared value of a^j is given by 



2sl 



dfl 



{alfp{a'^,)da^„=sl 



Using Eq. (l2Tl i resulting from equipartition of the internal 
bending energy, the standard deviation is given by 



-q-' (for«>l). 



In the absence of an externally applied force, a 2D semiflexi- 
ble filament characterized by the tangent angle 9 



0(5,K}) = ,/— £fl«COsM, (C-1) 

is in a configurational state denoted by {aJJ}. The probability 
of a filament being in a state between {a^} and djaj]} is 



p(K})dK} = 2^'^p 



knT 



Y{K. (C-2) 



F7>1 



where djaj)} = nn>i daj], Z is the partition function given by 



exp 



n>\ 



and is the internal bending energy of the filament given by 
Eq. ( I20I 1. After substituting the Fourier series for Q into 
and by using sin w 0, one finds 



_ Lp 



(C-4) 



;i>l 



This result can be substituted into Eq. (|C-3t yielding the par- 
tition function 



Z = Y[ /exp 



I7>1 



'P 2/ Q\2 



„>l \Lpqi, 



(C-5) 



The ensemble-average of a quantity AdajJ)}) is defined by 



(A) 



p(K)})A({fl«)})nd«° (C-6) 



I7>1 



The end-to-end distribution function can be calculated by the 
following ensemble average^^ 



' n>l 



(C-7) 



where 8{ro) is the Dirac delta function, and Lq — j L«>i ('^n)^ 
is the end-to-end distance of the chain in the absence of an 
externally applied force (/^ = 0), according to Eq. ( fT6] l. Next, 
we can combine Eqs. ( |C-2|C-4|C-7t and make use of the fol- 
lowing property of the delta function 



5(ro) = ^ / exp[/zro]dz, 



resulting in the following integral solution for $#(ro) 
^^(ro) = i- / dzexp [iziro - Lc)] ^ 



1/2 



(C-8) 
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Next, we rewrite the product inside the integral expression as 

1 



n 



/ «2 \ 1/2 



„>l \Lpql-iz 
in which \j/ is defined as 



: exp 



(C-9) 



yf{z) = 1^ In 



„>i 

The derivative dxir/dz can be written as 



dz 



1 



1 



^iLpqf;-z 2z 



1 




cot 




and integrated to yield (with ^/(O) — 0) 
1/(2) = In 



zLl/Lp 



(C-10) 



sin J zL^/L-p 



Equations (IC-9IC-10b are subsituted into Eq. dC-SI ), and by 
substitution of variables accoring to z = iLpCO^ we eventually 
arrive to 



'^{i-q) = ^ / doexp [-~(O^Lp{ro-Lc)] 



X 2a ( ] 
Vsinh(toLc)/ 



1/2 



(C-11) 



Next, we can use the binomial expansion^ to expand the in- 
verse square root of sinh(a)Lc) 



sinh((0Lc)"2 =%/2^ 2 (_i)*exp [-(2fc+ i)toLc] 



A'>0 

The distribution function then becomes 
LpV2L^ 



X exp 



-(O^Lp{ro ~Lc)- {2k + -)(0Lc 



(C-12) 



We can rewrite Eq. (IC-12b using so-called paraboUc cylinder 
functions Dy (z), defined by^ 



/V 7 

By (x) = exp — / i exp —xs + —s ds. (C-13) 



1 


V 




1 2" 


, — exp 
V27r/ 


4 


exp 





Using the following substitution of variables, 

, -, 2(k+l /4)Lc 

s = ft)V2Lp(Lc - ro), andx= =, 

y^2Lp(Lc-ro) 



the distribution function is rewritten as (with v = 3 /2) 

1 



^('"0) = — ^ L 

k>Q 



(-1) T 



X exp 



2Lp(Lc - ro) 



D 



3/2 



2(fc+|)Lc 
^2Lp(Lc - ro) 



(C-14) 



Finally, we can rewrite the binomial coefficient in terms of the 
Gamma function: 



r(fc+i)r(i-fc)' 



with r(A;+ 1) = By making use of the following identity 
of the Gamma function 



r « 



1 



(2n-l)!! j./l 
2" V2 



substituting n~—k and the identity 

(-1)' 



(-2A:-1)!! = 



{2k-\)\V 



we arrive at the the following expression for the binomial co- 
efficient 

-h\_ {2k-\)\\ , 
k )^ 2'kl ^ ' ■ 

The distribution function is therefore written as 

2LpVIc (2;t-l)!! 1 



X exp 



<r>0 



2kk\ [2Lp(Lc-ro)]5/4 



2Lp(Lc-ro) 



D 



3/2 



2{k+\)Lc 
yj2Lp(Lc-ro) 



(C-15) 



Appendix D 

We consider a 'network' consisting of parallel filament strings 
that are initially normal to the plane of shear, separated by a 
distance H [see dashed lines in sketch(a)]. The strings have 
a stretching stiffness /i and connect the top and bottom plates 
separated by a distance W . Next, the network is subjected to 
a shear of strain Y by displacing the top plate by a distance 
u. Strings rotate by an angle )3 = arctan(r) and undergo a 
stretch, leading to a certain shear stress t. The instantenous 
stiffness G can be found by increasing the top-plate displace- 
ment from M to M + Am [see sketch(b)] corresponding to a in- 
crease AF in strain and resulting in an increase At in stress. 
Using a unit out-of -plane dimension, the shear stiffness is de- 
fined as 



At _W M"^ 
~ AT ^ 77 aI7' 



(D-1) 
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(b) 



1/1/ W 



u Au 



L//L+AL 




Inserting Eq. (ID-3b into Eq. (ID- lb and using u = Wtanj3, we 
arrive at 



G = — sin^ j3 cos j3 . 
H 



(D-4) 



in which APx is the increase in string's force (AP) projected 
on the shearing direction, i.e., APx ~ APsinp. The linear re- 
sponse of each string. 



Finally, since the distance between strings is inversely pro- 
portional to the network's line density, i.e., H °^ p^^, and by 
employing j8 = arctan(r) we obtain at the following scaling 
relation for the stiffness 



AL 



(D-2) 



in terms of its elongation AL, along with the geometrical rela- 
tions AL = Am sin j5 and L = u/ sin j3 , yields 



APx = IJ. — sm^jS. 
u 



(D-3) 



PPL 



(i+r2)3/2 



(D-5) 
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